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We use a perturbative renormalization group approach with short-range continuum model inter- 
actions to analyze the competition between isotropic gapped and anisotropic gapless ordered states 
in bilayer graphene, commenting specifically on the role of exchange and on the importance of spin 
and valley flavor degeneracy. By comparing the divergences of the corresponding susceptibilities, we 
conclude that this approach predicts gapped states for flavor numbers N = 1, 2, 4. We also comment 
briefly on the related gapped states expected in chiral (ABC) trilayer graphene. 
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I. INTRODUCTION 

In chirally stacked few layer graphene systems 1 - - — two 
sublattice sites, one in the top layer and one in the bot- 
tom layer, are disconnected from the near- neighbor in- 
terlayer hopping network and partially isolated. This ge- 
ometric feature leads to a gapless semiconductor with 
weakly dispersive conduction and valence band edges, 
and therefore to an enhanced role for electron-electron 
interactions^. The low-energy sublattice degree-of- 
freedom in chirally stacked N layer graphene can be 
viewed as a S = 1/2 pscudospin, with t-pseudospins 
corresponding to top layer states and 4,-pseudospins 
to bottom layer states. In this language the conduc- 
tion and valence band eigenstates at wavevector q = 
q (cos 4> q , sin <p q ) have x—y plane pseudospin orientations, 
i.e. they have equal weight in top and bottom layers, 
and azimuthal angles N<j> q + it and N<p q respectively. A 
number of years ago we predicted^ on the basis of an un- 
restricted Hartree-Fock mean-field-theory that these cir- 
cumstances would lead to ground states with pseudospin 
order, and in particular to a state in which the pseu- 
dospins rotate out of the x — y plane by a fc-dependent 
polar angle toward either the north or the south pole, 
breaking layer-inversion symmetry 7 - and inducing a gap£ 
in the charged excitation spectrum. Recent bilayer ex- 
perimental worker— has yielded clear evidence for broken 
symmetry states, although a consensus on the character 
of their order has not yet been achieved. At the same 
time a substantial body of theoretical work^&i£~— has 
used a variety of different approaches to study bilayer 
graphene. The general consensus of this work is that the 
ground state is a pseudospin ferromagnet. It has become 
clear however that because of quantum fluctuations that 
are not captured by mean-field theory, the gapped state 
competes closely with a gapless anisotropic state in which 
the net pseudospin alignment is in the x — y plane, not 
the z direction. This paper addresses that competition. 

As we will explain in more detail later, some features 
of the bilayer graphene system present awkward obsta- 
cles to theory. First among these is the property that the 
conduction and valence bands arc weakly dispersive only 



over the small part of the Brillouin zone where inter-layer 
tunneling plays an essential role. Elsewhere in the Bril- 
louin zone interaction effects arc expected to be similar 
to those in single-layer graphene, producing strong quasi- 
particle velocity enhancements ^Sr— and spectral weight 
rearrangements ; 41 ' 42 but not order. The small number 
(~ 10~ 4 ) of 7r-electrons per carbon atom in the strongly 
correlated region of momentum space is an obstacle for 
theories of the ordered state that are based primarily on a 
lattice model of graphene. The narrow momentum-space 
distribution of the most important virtual states also cre- 
ates difficulties for continuum models because it makes 
the physics sensitive to the long-range of the Coulomb 
interaction. In this paper we report on calculations in 
which lattice effects are completely ignored and the long- 
range of the Coulomb interaction is not treated explicitly 
In our view the continuum approach we use is strongly 
motivated by the low-density of strongly correlated elec- 
trons. Our use of short-range interactions can be crudely 
justified by appealing to screening considerations, and 
arguing that the momentum-independent interaction pa- 
rameters we use represent an average over the relevant 
portion of momentum space. Although not fully rig- 
orous in this sense, we believe that the conclusions we 
reach by following this approach systematically are nev- 
ertheless valuable and consistent with recent experimen- 
tal work^iM 

This paper elaborates on ideas presented in a series 
of closely related earlier papers^ ' 17 ' 19 ' 22 ' 31 in which the 
ground state is estimated on the basis of perturbative 
renormalization group (PRG) calculations. These calcu- 
lations differ among themselves most essentially in the 
approximate ways they account for lattice scale physics 
and long-range interactions. In Section II we briefly de- 
scribe the model we study and restate some identities 
and observations that are useful in developing its per- 
turbation theory. Section III presents new results for the 
PRG interaction flows of this model, emphasizing the role 
of exchange interactions between electrons with the same 
spin and valley states and commenting on the flavor num- 
ber (N) dependence of the interaction parameter flows. 
In Section IV wc explain how the PRG pscudospin sus- 
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ccptibilitics, whose divergences suggest the character of 
the low-temperature order, depend on the renormalized 
interaction strengths. Finally in Section V we present 
our conclusions and discuss the relationship of our work 
to experiments and to other theoretical work on this sys- 
tem. 



II. THEORETICAL PRELIMINARIES 

Bilaycr graphcne (BLG) is described approximately at 
low-energies by the 3 = 2 version of the chiral band 
Hamiltonianii — 

Hj = J-i c Ui I cos{J(j) q )al p +sm{J(j) q )a v a 



(i) 

In Eq. ([1]) <r are Pauli matrices which act on the (Greek) 
layer labels and i = 1, . . . , N is a spin- valley flavor la- 



bel. We have used the notation cos< 



z q*/q and 



sin</>g = q y /q. Here r z = ±1 labels K and K' valleys 
located at inequivalent Brillouin zone (BZ) corners and 
q = (q x ,q y ) is wavevector measured from the BZ cor- 
ner. It is convenient to perform a rotation in BLG pseu- 
dospin spaced in the K'(t z = —1) valley to eliminate 
the valley dependence of the band Hamiltonian, in order 
to make our discussion more concise in the following sec- 
tions. The BLG chiral band Hamiltonian applies^ at en- 
ergies smaller than the interlayer hopping scale 71 ~ 400 
meV and larger than the trigonal warping scale ~ 1 
mcV. We discuss the important role of trigonal warp- 
ing in the band-structure later. The physics discussed 
in this paper applies to bilayers when J = 2 and partly 
generalizes to J-layer chirally (ABC) stacked graphene 
multilayers ) i 43 ' ' although the chiral band Hamilto- 
nian Hj applies over narrower ranges of energy for larger 
J. We will focus our attention on the J = 2 case here- 
after. 

H2 can be viewed as specifying a momentum q- 
dependent effective magnetic held that acts on the bilayer 
layer-pseudospin degrce-of- freedom: 



n 2 



B er 



(2) 



where \B\ = £ q = hiv q = fi 2 q 2 /2m* with to* = ji/(2v 2 ), 
and the orientation angle of B is 2<fi q . The corresponding 
Matsubara Green's function is 

G q (iuJn) = [iuJ n -'H/K\^ 1 = Gqs,{iun)+G q t{iu n ) n a (3) 
where 



Sqs,t(iw n ) 



1 



1 ± 1 



2 V iu B — LO q icj n + 0J q 



(4) 



and n = — (cos 2<j) q , sm2<f> qi 0). Note that G s {—iu n ) = 
~G s (iu> n ) whereas Gt{— iw n ) = Gt(iw n ). When expressed 
explicitly as a 2 x 2 matrix, 



G s (iuj n ) ~Gt(iuJn)e- 2i ^ 
-Gt{i^n)e 2i ^ G s (iuj n ) 



(5) 



The off-diagonal triplet component captures processes 
in which electrons propagate between layers; its 
momentum-orientation dependent phase factor plays an 
essential role in the one-loop PRG calculations described 
in the next section. The frequency sums which appear in 
loop diagrams are readily evaluated: 

1 X-C 2 tanh(/^ q /2) 1 
w ^ G qs , t M = T ^ T^- , 

w^^K)^W^ ' ( 6 ) 

These two results express the property that fluctuations 
involve interband transitions. In the PRG each loop di- 
agram is multiplied by appropriate interaction constants 
(discussed below) and then integrated over momentum 
labels near the model's flowing cutoff A: 



A/s<g<A 



d 2 q 



tanh(ggg/2) 
4^ 



1 

T-s-0 2 



^0 ln(s). (7) 



where vq = to* /2ttH 2 is the BLG density-of-states per fla- 
vor. The second frequency sum in Eq.® vanishes in the 
long-wavelength static limit because it has contributions 
from intra-band particle-hole excitations only. Because 
uj q oc q 2 rather than q, this integral diverges logarith- 
mically when the high-energy cut-off is scaled down by a 
factor of s. The PRG calculation is therefore very similar 
to the corresponding calculation for a IDES. This rather 
surprising property of BLG is directly related to its un- 
usual band structure which yields Fermi points rather 
than Fermi lines in neutral systems, and quadratic rather 
than linear dispersion^ 



III. RG FLOW EQUATIONS 

A. Number of Physical Parameters 

In systems with short-range interactions, finite mo- 
mentum (gradient) corrections are irrelevant in the RG 
flows. We therefore ignore the wave- vector transfer de- 
pendence of scattering processes starting already at the 
J = 2 chiral model's ultraviolet cut-off energy ~ 71. In 
the low-energy continuum model of BLG electrons carry 
spin, and both layer and valley pseudospin labels. In a 
scattering event, both the two incoming and two outgo- 
ing particles can therefore have one of eight labels and 
the general scattering function therefore has 8 4 distinct 
values, even when gradient corrections are ignored. The 
number of distinct coupling constants in the RG flow 
equations is much smaller, however, because many val- 
ues are zero and others are related to each other by sym- 
metry. One simplification is that interactions conserve 
spin and, in the continuum model, both layer and val- 
ley pseudospin at each vertex. The role of spin in our 
model is therefore equivalent to the role of valley. We 
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do however choose to allow the interaction strength be- 
tween electrons in different (D) layers to be weaker than 
that between electrons in the same (S) layer so that the 
interactions are layer-pseudospin dependent. 

The internal loops in the PRG calculation contain two 
fcrmion propagator (Green's function) lines. Because we 
seek a momentum-independent effective interaction we 
simplify the calculation by evaluating the coarse-graining 
correction to the effective interaction for the case where 
all incoming and outgoing electrons are at the Dirac 
point. (As we comment in the discussion section, this 
approximation is made less innocent by the layer pseu- 
dospin structure of the band eigenstates.) These prop- 
agators conserve both spin and valley pseudospin, but 
as we have seen above, not the layer pseudospin. From 
Eq. ([5]) we see that a phase factor— e ±2l(?iq appears when 
a propagator transfers electrons between layers. Unless 
these layer transfers enter an equal number of times in 
each direction, the integrand in a Feynman diagram will 
contain a net phase factor related to chirality and van- 
ish upon angular integration of the virtual momentum. 
The total layer number (the z component of layer pseu- 
dospin) is therefore also conserved^ This property com- 
bined with flavor invariance limits the number of inde- 
pendent interaction parameters to three, two for inter- 
action processes involving electrons in the same layer or 
different layers that conserve layer index, and one for in- 
teraction processes between electrons in different layers 
accompanied which interchange layer indices. We refer 
to the three distinct processes, illustrated in Fig. [TJ as 
Ts, Td and Tx respectively. 



t T T T ! I 




t T 1 I ! I 



FIG. 1. (Color online) Electron-electron scattering processes 
for a system with a single pseudospin- 1/2 degree of freedom. 

If flavor invariance is broken, either due to microscopic 
lattice effects that violate SU (4) invariance or due to bro- 
ken symmetries, the number of interaction parameters 
increases. For example, if the bare interactions are spin 
and valley dependent in an arbitrary way, the number of 
independent antisymmetrized interaction parameters in- 
creases from three to ten£ Tssd^sds, Tsdd, r DSS ,r D sD, 
Tdds, r DDD , Txsd, Txds and Txdd- In this notation the 
first label follows the convention explained above for the 
layer pseudospin dependence of the interaction, whereas 
the second and third labels refer to valley and spin respec- 
tively, and distinguish the interaction between particles 
with the same (S) or different (D) spin or valley. The 
number of entries in this list of interactions has been re- 
duced by two by appealing to Fermi statistics for the case 



of SS interactions, those between particles with the same 
spin and valley labels. In that case interchanging the la- 
bels of outgoing electrons yields an indistinguishable and 
canceling amplitude for particles in the same layer. Al- 
though we can choose to retain Tsss in the RG flows, it 
cannot contribute to any physical observable. Similarly, 
since Tdss and Txss are interchanged by reversing out- 
going labels they differ only by a sign and it is sufficient 
to retain one parameter. The number of independent in- 
teracting parameters in the analysis by Vafck 2 ^ and by 
Lemonik et aZ.)24 who impose the symmetry of the un- 
derlying BLG crystal but do not explicitly account for 
antisymmetrization is nine. For SU (4) invariant models 
the number of parameters is three, but antisymmetriza- 
tion allows us to retain only one coupling constant in the 
single flavor (N = 1) case, and the physical implications 
of the flowing interactions are more apparent if we do 
so.— A convenient choice^ is to retain Tdss- We take the 
view that the underlying lattice 2 ^ is unlikely to play an 
important role in determining whether the broken sym- 
metry state is gapped^ or gapless , 17 ' 22 though it is key 
to selecting 2 ^ between distinct gapped states^ that are 
equivalent in an SU (4) invariant model. Therefore, the 
number of interactions parameters that appear in the cal- 
culation described below is three in the general case with 
unbroken flavor invariance and one for the special case of 
N = 1. 



B. RG Flow for N = 1 

In this subsection we explore the similarities and differ- 
ences between BLG and IDES's by temporarily neglect- 
ing the spin and valley degrees-of-freedom. We choose 
the single effective interaction parameter to be Td- A 
PRG analysis determines how the bare interaction at the 
7i scale, Vb, is renormalized by integrating out the high 
energy fermion degrees of freedom. At one loop level 
the coarse-graining contributions to the effective interac- 
tion are described^ 7 - by the three higher order diagrams 
labeled ZS, ZS', and BCS. The internal loops in these 
diagrams are summed over the high-energy labels. The 
main merit of the PRG is that it treats all virtual pro- 
cesses on an equal footing, 

It is well known that in a IDES the ZS loop vanishes 
while the ZS' and BCS diagrams cancel, implying that 
the interaction strength does not flow to large values 
and hence that neither the CDW repulsive interaction 
nor the BCS attractive interaction instabilities predicted 
by mean-field theories survive the quantum fluctuations 
they neglect. The key features of BLG order physics 
can be understood in terms of two properties of these 
one-loop diagrams; (i) the particle-particle (BCS) and 
particle-hole (ZS, ZS') loops have the same logarithmic 
divergences as in the IDES case in spite of the larger 
space dimension and (ii) the ZS loop, which vanishes in 
the IDES case, is finite in the BLG case and the BCS 
loop vanishes instead. Both of these changes are due to 



the layer pseudospin triplet contribution to the single- 
particle Green's function as we explain below. The net 
result is that interactions flow to strong coupling even 
more strongly than in the mean-field approximation. 

The key step in one-loop PRG calculations is identify- 
ing the coupling factors attached to the loop in each di- 
agram. Since only opposite layer interactions need to be 
retained for N = 1 models, all scattering functions have 
two incoming particles with opposite layer labels and two 
outgoing particles with opposite layer labels, much like 
the IDES case4£ The external legs in the scattering func- 
tion Feynman diagrams are labeled by the layer index 
(T = top layer and B = bottom layer). The correspond- 
ing labels for the IDES cas o 47 i 48 are the chirality, R for 
right-going and L for left-going. We refer to these as 
the single-particle labels below when a comment refers 
to both IDES and BLG cases. 



(a) 



(b) 





FIG. 2. (Color online) (a) ZS, (b) ZS' loop corrections in the 
one-loop PRG calculation. 

As shown in Fig.[2ja), at the upper vertex of the ZS di- 
agram the incoming and the outgoing T particles induce 
a B particle-hole pair in the loop, while the incoming 
and outgoing B particles at the lower vertex induce a T 
particle-hole pair. The ZS contribution is absent in the 
IDES cas e 47 i 48 because propagation is always diagonal in 
interaction labels. However, this correction survives for 
BLG because the single-particle Green's function has a 
triplet contribution (see Eq. (|6])) which is off-diagonal in 
layer index. Here we find 



r 



zs 



0/(^E^-n) = ir^ o ln( S ).(8) 



The ZS' loop shown in Fig. EJb) corresponds to re- 
peated interaction between a T particle and a B hole. 
This is the channel responsible for the IDES mean-field 
CDW instability^ 7 - in which coherence is established be- 
tween R and L particles^ 8 - In both IDES and BLG 
cases it has the effect of enhancing repulsive interac- 
tions. Its evaluations for the two cases correspond quite 
closely, because this loop diagram involves only particle- 
propagation that is diagonal in interaction label. We find 
that 
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FIG. 3. (Color online) BCS (particle-particle) loop correction 
for singlet propagation in the one-loop PRG calculation. 



The BCS channel corresponds to repeated interaction 
between the two incoming particles. In the IDES case 
the contribution from this loop (see Fig. [3]) cancels the 
ZS' contribution 4 - 7 -, leading to marginal interactions and 
Luttinger liquid behavior. This same kind of BCS cor- 
rection for BLG reads 



r 



BCSi 



1 T 2 D f d 2 q ^ . 

i(-r D ) 2 f d 2 q ^ 

J (2^)2 Z^G s {q^^)G s {-q,-i^n) 

(10) 



2 f3h 2 

1 



- ^F 2 ) v ln(s) 




FIG. 4. (Color online) BCS (particle-particle) loop correction 
for triplet propagation in the one-loop PRG calculation. 

In the BLG case there is an additional BCS loop con- 
tribution (see Fig. [4j> in which the incoming T and B 
particles both change layer labels before the second in- 
teraction. This contribution is possible because of the 
triplet pseudospin propagation and, in light of Eq. (|6|), 
gives a BCS contribution that has the opposite sign com- 
pared to the normal contribution: 



F 



BCS 2 



ir D 

2" 



r D ) f d 2 q 



1 



j3h 2 

1 (-r D )r D 

2 f3h 2 
r 2 ) v ln(s) . 



(27T) 2 

d 2 q 



Gt (<?, iu n )Gt (-9, -iu) a ) 

^ Gt (q, ^ n )5t (-<?, -iWn) 

(11) 
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TABLE I. Summary of contributions from the three one-loop 
diagrams to Td in units of v in IDES and BLG cases. 



tion for BLG: 



diagrams 


1DEG 


BLG 


ZS 





| T 2 ln( fl ) 


ZS' 


u 2 ln(s) 


- Vl hits) 


BCS 


-u 2 ln(s) 





Mean Field 
Quantum Fluctuations 
Full One-Loop 


u 2 ln(s) 
-u 2 ln(s) 



ir 2 ln( fl ) 
I rl ln(s) 
T 2 ) ln(s) 



It follows that the BCS loop contribution is absent in 
the BLG case because 



t^BCS t^BCSi | t^BCS2 n 

1 D — 1 D +t D — U ' 



(12) 



Therefore, at one-loop level, the renormalization of inter- 
layer interaction is 



^one-loop _ r ZS , -pZS' , pBCS _ y.2 



D 



rg> + rg> +rg^=rl,^in(3). (13) 



These results and the comparison with IDES's are sum- 
marized in Table U and imply the following RG flow equa- 



i>o rfln(s) 



- F 2 



(14) 



If we follow the flows from the model's microscopic cutoff 
we obtain 



r n = 



Vd 



1 - V D i / oln(s) 



(15) 



where Vn is the bare interaction. Tn diverges when 
Vy) > l/ln(s). This instability criterion is similar to 
the Stoner criterion for ferromagnetism. 

The interaction correction to the layer pscudospin re- 
sponse function Xzz is obtained by closing the scattering 
function with a <r z vertex at top and bottom. The a z 
operator measures the charge difference between T and 
B layers. Because it is an effective single-particle theory, 
fermion mean-field theory corresponds to response func- 
tion diagrams with at most a single particle-hole pair. 
It follows that mean-field theory in BLG is equivalent 
to a single-loop PRG calculation in which the BCS and 
ZS' channels are neglected and only the ZS channel is 
retained. In mean-field theory, ideal BLG has an insta- 
bility to a state in which charge is spontaneously trans- 
ferred between the layers which is signaled by the di- 
vergence of % zz . (xxx also diverges in mean- field theory 
but less strongly.) The PRG analysis demonstrates that 
the mean-field theory instability is enhanced by the re- 
inforcing influence of the ZS' channel contribution. We 
discuss pseudospin response functions which indicate the 
character of the broken symmetry state further for the 
physically relevant N > 1 case in Section IV. 



RG Flows for N = 1, 2, 4 



Assumming SU (4) invariance, the one-loop flow equations for N > 1 are derived in the same way as for the 
N = 1 case, except that we need to keep track of three interaction parameters, Tg (Tg = Tsds = Tssd = Lsdd), 
Td (I'd = Tods = Tdsd = T D dd) and T x (T x = Txds = Txsd = T X dd)- A tedious but elementary book-keeping 
exercise yields the following results: 



s 



r§ - (r D - F x )(r D - r s ) - i— -(r D - r s f + ^(r x - r s ) 2 , 

dT B _ 1 2 ,„ „ w „ „ , N-2,„ _ 1 



r D + (r D - r x )(r D - r s ) + -^^(r D - r s ) 2 - -(r D + r x )' 



vq rfln(s) 2 2 2 

dVv N-2 9 1, N , I, x9 

= (r D - r x )r x - — — r x - -(r x - r s ) 2 - -(r x + r D ) 2 . (16) 



uodln(s) y " ~' A 2 2 °' 2 

I 

This equation is the SU(4) invariant version of the more in the combination Td — T x . Defining Tds = Ld — T x , 
general ten parameter flow equations derived in Ref. H|. it follows from Eqs. pi)]) that 



The RG flow equations derived previously in Ref. 1171 differ 

only in notation. The N = 1 case is recovered by not- dTns 2 N — 1 . 2 N — 1 2 

ing that T g has no physical effect for interactions among z/ dln(s) ~~ ^ DS ^ 2 _ ^ S ^ 2 
particles of the same flavor and Td and T x always enter 



X ) 

(17) 
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yielding a one-parameter flow equation for the N = I 
case. Note that this Trjs also represents rrjss in the 
N = A case.- 

The only fixed point of these RG flow equations is the 
non-interacting one. The character of the broken symme- 
try state can be estimated by following the RG flows to 
strong interactions, or more reliably by following the flow 
equations for appropriate susceptibilities as discussed in 
the next section. In order to represent the property that 
bare same layer interactions should be slightly stronger 
than the bare different layer interactions and that the 
bare layer interchange interaction vanishes, we start the 
RG flow integrations from 



V S >V D , Vx = 0. 



(18) 



When the flows start from an interaction with these prop- 
erties, the coupling constants diverge simultaneously at a 
finite value of ln(s), and for N = 2, 4 satisfy the following 
properties at the divergence point: 



r D 
£s 
r x 

r x 



-0.5425 (for N = 2) , 
-0.2624 (for N = 4) . 



(19) 



The value of the first ratio in Eq. (|T9| follows from the 
RG flow equation for Ts + Td which satisfies 



rf(r s + r D ) 

vq dln(s) 



-r x (r s + r D ). 



(20) 



The rate of growth of Ts + Td is proportional to |r x |, 
whereas Ts and |T X | grow like |T X | 2 . Equivalent re- 
sults in a different notation, using 2go,z = Ts ± Yd and 
2g± = r x , have been derived previously by Vafek and 
Yangi^i These results do not depend on the bare coupling 
constants as long as the physical conditions specified by 
Eq.(fI5]l are satisfied. Note that the RG flows for N = 4 
are still qualitatively different from those of the N —> oo 
limit. 

Typical RG flows, obtained by integrating Eq. dTBl nu- 
merically, are plotted in Fig. [5] We note that, when com- 
bined with the mean-field-theory of the gapped state^^ 
the value isqVs = 0.25 yields values of the spontaneous 
gap, the critical temperature, and the carrier density at 
which order is absent that are consistent with some re- 
cent experimental observations^^. We therefore choose 
to use the bare interaction parameters vqVs = 0.25 and 
voVd = 0.22 to illustrate typical RG flows. We find 
that the interaction parameters flow away from the non- 
interacting fixed point and diverge at a finite value of 
s. The instability criterion implied by the one- loop RG 
calculation is vqV$ ~ 0.7/ln(s) for N = 4, i..e the in- 
stability tendency is enhanced compared to the N = 1 
case. (Although the Tg, Trj and Tx interaction parame- 
ters do diverge at a similar value of ln(s) in the N = 1 



case, they have no physical effect because they enter ob- 
servables in the combination Trjg = Trj — T X -) Because 
of the quadratic band dispersion, the ratio of the inter- 
layer coupling scale to the gap (71 / A) can be associated^ 
with the square of the momentum scaling factor s at the 
phase transition. This yields a value of A on the order of 
1 meV, close to the scale of the gap values seen in some 
experiments! 9 ' 10 ' 12 " — 



(a) N=1 




(b) N=2 
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(c) N=4 
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FIG. 5. (Color online) RG flows of the interaction couplings 
for (a) N = 1, (b) N = 2 and (c) N = 4. In the N = 1 
case the dashed lines are unphysical since interactions always 
enter observables in the combination Tds = Td — Tx- 



IV. SUSCEPTIBILITIES 



More insight into the likely nature of the broken sym- 
metry state which occurs in BLG can be obtained by 
using the PRG calculations to estimate the long wave- 
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FIG. 6. (Color online) Feynman diagrams for the pseu- 
dospin susceptibilities, (a) The non-interacting susceptibil- 
ity and (b) the interaction correction to the susceptibility. 
T,B — x,y, z label layer-pseudospin components. The spin- 
valley flavor labels are not explicitly indicated. 



length static limit of the pseudospin susceptibilities: 

. 

X t 7 ',b 7 {q, iw)= dre^ T (T T S T y (q, t)Sb 7 (-9, 0)) , 
Jo 

(21) 

where T, B = x, y, z label layer-pseudospin components 
(instead of top and bottom layers) and 7', 7 label either 
a unit matrix (7 = 0) or a generator of rotations (7 = 
1, . . . , 15) in the 4-component valley-spin flavor space: 

S^(q,T)=J2 c l+ q , a i>°S A ?'i c W- ( 22 ) 

k 

A physically transparent choice^ for the 517(4) gener- 
ators is to set AJ, i to a Pauli matrix which acts on the 
spin dcgrcc-of- freedom for 7 = 1,2,3, to a Pauli matrix 
which acts on the valley degree-of-frcedom for 7 = 4, 5, 6, 



and to a product of spin and valley Pauli matrices for 
7 = 7, . . . , 15. It follows from SU(4) invariance^ that the 
normal state susceptibilities in this case are oc <5 7 < i7 and 
equal for all values of 7 > 1 . Below we refer to the 7 = 
susceptibility (x c ) as the flavor charge susceptibility and 
to the 7^0 susceptibility (x s ) as the flavor spin suscepti- 
bility. The conservation of particle-number in each layer 
at long wavelengths, explained previously, implies that 
Xt 7 ,b 7 oc St,b and that Xx 7 ,x 7 = Xy 7 ,y 7 - Divergences 
in z signal ordered states in which the charge density 
is transferred between layers in the same sense for all 
flavors, while divergences in xl z signal flavor dependent 
layer-inversion symmetry breaking with no overall charge 
transfer. Divergences in Xxx f° r T = x,y, on the other 
hand imply broken symmetry states in which the phase 
relationship between layers is altered and rotational sym- 
metry within the layers is broken. Broken layer-inversion 
symmetry states have an energy gap for charged excita- 
tions and large momentum space Berry curvatures that 
lead to large anomalous Hall effect contributions from in- 
dividual flavors^ These states will be referred to below 
as spontaneous quantum Hall states^ The broken rota- 
tional symmetry state will be referred to below as the 
ncmatic state. 

Provided that the system has a single continuous phase 
transition, the ordered state character can in principle be 
determined by identifying the pseudospin susceptibility 
which diverges first as temperature is lowered. Although 
we cannot make that determination definitively from the 
RG interaction strength flows, we can obtain some guid- 
ance by examining how the interaction contributions to 
the susceptibilities grow under the RG flows. The pseu- 
dospin response functions of interest can be related to 
the renormalized interactions as summarized in the Feyn- 
man diagrams of Fig. [S] Fig. represents the non- 
interacting pseudospin susceptibility x° and Fig. |6fb) 
represents the interaction correction x' nt - I n the long 
wavelength static limit the explicit expressions for these 
susceptibilities arc: 



I 

X% = ~ [ (2 J 2 ^ h2) Y, J2 ^ea T ('7,«^)^ ) (3T ^ T(3B (g,^) ( T2 ) QB , (23) 

^ ' ^ ' u q t /3 t «b/3b 

Xtb = ~ [(2 n )4(3fj?)2 ^ ^ ^aT^i^iV^^^ (24) 



Note that in these equations all propagators are diagonal 
in the implicit spin- valley flavor labels. 



In the noninteracting case, the inter-flavor (D) suscep- 
tibilities vanish because the single-particle Hamiltonian 
is diagonal in flavor and the intra-flavor (S) susceptibil- 
ities are positive because of the band state pseudospin 



structure: 

OD OD _ yOD q 

Azz a xx ^yy ' 

XL S = 2x^ = 2 X °^ 2^ ln( S ). (25) 

These pseudospin susceptibilities capture the contribu- 
tion of momentum q = vertical interband quantum 
fluctuations involving states with energies measured from 
the Dirac point between 71 and 71 /s 2 . The factor of two 
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difference between Xzf and Xxx demonstrates that the 
band state is more susceptible to a gap opening pertur- 
bation than to a nematic perturbation. This property 
can be understood 2 ^ in terms of the pseudospin orien- 
tations of the band states; in particular a z pseudospin 
effective field is perpendicular to the valence band pseu- 
dospin for all momentum orientations, so that I direction 
fields always yield a strong response. On the other hand 
x — y-planc pseudospin fields are not in general perpen- 
dicular to the valence band pseudospin and consequently 
produce a weaker response averaged over momentum ori- 
entations. The larger response to a z pseudospin effective 
field is related to the well-known property of BLG that 
a potential difference between the top and bottom layers 
lead to an energy gap at the Dirac point. The corre- 
sponding charge susceptibility (density-density channel), 
obtained by replacing the pseudospin-ficld Pauli matrix 
by an identify matrix, vanishes because conduction and 
valence band states are orthogonal. 

When the interactions flow to strong values, the sus- 
ceptibilities are dominated by their interaction contribu- 
tions. When these are included we find that 

xZ = 2(^ in( s )) 2 (r D -r s ), 
x^-^Ki^frx, (26) 

X S ZZ = 2^o ln(s) + 2(^o ln(s)) 2 (r D - I*) , 

XL = W ln( S ) + ifo ln( S )) 2 (r D - T X ) . (27) 

Note that the large susceptibility in Xxx is attribute to 
the layer interchange processes Tx while the strong di- 
vergences in xtz xx are due to intra-flavor exchange inter- 
actions, recalling Ld — Lx = Tds in Eq. (f27j) . 

For the N = 1 case there is no different-flavor response 
and we therefore find that near the instability 

Xzz = 4Xxx = 4 Xyy , (28) 

suggesting that a gapped spontaneous quantum Hall 
state is most likely. 

For the physical N = 4 case we can relate the flavor 
charge and flavor spin response functions to x S and x D by 
examining the particular case of responses to perturba- 
tions that are diagonal in flavor index, i. e. perturbations 
proportional to s z , t z or s z r z . We find that near the 
divergence point 

Xzz J Xzz ^ Xzz J (29) 
Xzz ^ Xzz ^ Xzz ' 

(30) 

Xxx have the same relations to Xxx and Xxx- Combining 
Eq.® and Eq.([27]) with the results from Eq.^ for 



N = 4wc find that 

X xx -> +1.7805 r s Kln( S )) 2 , 

Xzz ^ +2.4055 r s (^ln( S )) 2 , (31) 

X z c z ^ -1-5945 r s (^ln(s)) 2 , 

XL ~> -0.1250 Ts^olnW) 2 . (32) 

The property that all susceptibilities diverge together 
is an artifact of the single-loop RG calculation which 
becomes questionable when the scaled interactions are 
strong. Still, the fact that the xl z diverges most strongly 
suggests that the spin-channel gapped state is more likely 
than the charge channel nematic state. (Negative coef- 
ficients in Eq. ()32[) imply that one- loop interaction cor- 
rections tend to reduce the susceptibilities xL and Xxx 
rather than to increase them and therefore do not suggest 
instabilities.) 

In the N = 2 case, which accounts for only spin or 
valley but not both, XT7.B7' the 4 dimensional 7-matrices 
are replaced by 2-dimcnsional Pauli matrices. We find 
that 

X zz 2 Xzz ' 2 Xzz ' 

(33) 
(34) 

which equally applies to the Xxx channel. Combining the 
N = 2 results in Eq.(Q]§ with Eq.(gSJ) and Eq.^ we 
obtain 

XL^ -1-1567 rsKln(s)) 2 , 
XL^ +2.8433 r s (n,ln(«)) 2 , 
XL^ +0.6717 r s (^ ln( S )) 2 , 
xL-> -0.2500 r s (^ ln( S )) 2 . (35) 

In both N = 2,4 cases interactions favor breaking 
layer-inversion symmetry in the flavor spin channel and 
orientational symmetry in the flavor charge channel. The 
strongest divergence occurs for the flavor spin channel 
broken layer-inversion symmetry state.— The N = 1 sus- 
ceptibility and the flavor charge and spin channel suscep- 
tibilities for N = 2,4 are illustrated in Fig J7] and agree 
with the above analytic analysis. Because x zz i s always 
the most positvely divergent channel we conclude that the 
dominant many-body effect in BLG is spin-channel spon- 
taneous layer-inversion symmetry breaking yielding BLG 
states with spontaneous charge gaps. These results ver- 
ify our original theoretical predictions^ and are consis- 
tent with recent experimental observations in dual-gated 
ultra-clean suspended BLG devices^ ' 10 ' 13 i 14 

Finally, we point out that an interlayer potential 
breaks SU(4) invariance in the N — 4 case, leading to 
a range of stability for a state in which three flavors 
are polarized toward one layer and one flavor toward 
the other layer when an ordered state appears. For in- 
stance, a possible mass term in a mean-field theory with 
broken layer-inversion symmetry can be proportional to 
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(a) N=1 

10 1 1 1 




ln(s) 



FIG. 7. (Color online) RG flows of interacting suscepti- 
bilities for different ordering tendencies, (a) N = 1, (b) 
N = 2 and (c) N = 4. We use the dimensionless quantity 
^o(ln(s)) 2 (x T ) _1 in the figures. 



(1 + s z + t z - s z t z ) a z /2. This gapped state breaks par- 
ity inversion and time reversal symmetries, favored by 
small interlayer electric and magnetic fields^^ that break 
SU(4) symmetry. The difference between this state and 
other spontaneously broken symmetry states can be fur- 
ther appreciated as follows in our full model with ap- 
proximated 5?7(4) symmetry. This state reduces the full 
symmetry to SU{3) x SU(l) x U(l) with three Gold- 
stone modes, in contrast to the three possible states in 
flavor spin channel that break the symmetry down to 
5*7(2) x SU{2) x U(l) with four Goldstone modes, and 
there is no soft mode in the flavor charge channel. 



V. DISCUSSION 

In this paper we have attempted to shed light on 
the pseudospin ferromagnet ground state of BLG by 
carrying out a conventional perturbative rcnormaliza- 
tion group calculation for a model with 517(4) invari- 
ance in spin- valley space, and examining interaction cor- 
rections to spin and charge pseudospin susceptibilities. 
The PRG flows lead to strong layer interchange scat- 
tering processes which favor gapless nematic (x — y- 
plane pseudospin ferromagnet) ground states, compet- 
ing with the gapped (indirection pseudospin ferromag- 
net) states favored by intra-flavor exchange interactions 
and predicted^ by mean field theory. The role of ex- 
change is especially dominant in the single-flavor N = 1 
case in which the PRG flow equations lead to divergent 
interactions which vanish when projected onto the many- 
fcrmion Hilbcrt space. For N = 1 the interactions enter 
observable only in the combination Tds = Td — Tx that 
is uniquely allowed by the Pauli exclusion principle. 

We conclude that the most likely ground states are 
gapped spin-channel Ising pseudospin fcrromagncts, in 
agreement with somei2~— recent experiments, although 
the competition becomes closer when the number of fla- 
vors is larger, and quite close for the physically rele- 
vant N = 4 case. According to our estimates, broken 
symmetry states appear at energy scales slightly larger 
than those at which trigonal warping effects in the single- 
particle Hamiltonian^ start to have a large influence on 
the band structure, splitting the vorticity J = 2 vor- 
tex of the chiral model into four J = 1 Dirac points. 
The strength of our conclusions is limited by our neglect 
of the long-range character of Coulomb interactions and 
by the fact that our one-loop PRG calculations become 
unreliable once the effective interactions flow to strong 
values, i.e. once v$T > 1. 

The details of our calculations are conventional and 
in particular evaluate the loop corrections to scattering 
amplitudes for the case in which both incoming and out- 
going momenta are at the Dirac point. This simplifi- 
cation is normally justified by the irrelevance of gradi- 
ent corrections to the effective interactions. As we now 
explain, in the case of BLG, this approximation likely 
overestimates the strength of the Tx layer interchange 
interactions. Consider for example the ZS process which 
is non-zero because of the triplet layer off-diagonal contri- 
bution to the Green's functions, as we have approximated 
in Eq.([S]). When the incoming electrons have finite mo- 
menta k and p, the momentum arguments of the virtual 
states in the loop diagram are k-\-q and p+q. The triplet 
propogator phase factors which are integrated over q in 
the loop correction to Tx therefore becomes e ±2lS where 
9 = 4>k+q — 0p+q- Notice that this phase factor has a 
singular dependence on k and p at small q and that the 
approximation we have used previously, in which we re- 
place 8 by zero, is valid only for q 3> p,k. When the phase 
factors are taken into account the Tx scattering ampli- 
tude has an additional phase factor which will average 
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to a value smaller than one in physical obscrvablcs like 
the pseudospin susceptibilities. It therefore appears to 
us that our calculation likely overestimates the strength 
of the interactions which favor the nematic state. 

The gapped states have the interesting property 
that they support large momentum space Berry 
curvatures^i££ which can lead to interaction-induced 
quantum Hall effects. The regions of momentum space 
near both K and K' valleys flavor contribute ~ ±e 2 //i 
to the Hall conductivity for each spin. When all fla- 
vors make a contribution of the same sign, the total 
Hall conductivity is exactly ±4e 2 /h. This gapped bilaycr 
state is therefore a quantized anomalous Hall (QAH) 
insulato r 5 i 21 i 50 i 51 with chiral edge states. Because of its 
nontrivial overall orbital moment the QAH state is fa- 
vored by a perpendicular magnetic field. The state in 
which the sense of layer polarization is opposite for op- 
posite spins and for opposite valleys does not break time- 
reversal symmetry and has a mean-field interaction terms 
with an effective spin-orbit coupling^ that leads to a 
quantum spin Hall (QSH) state with zero charge Hall 
conductivity. The QSH state has helical edge states pro- 
tected by time-reversal symmetry. Its Z2 classification 
in an interacting iV-layer chiral graphene is given by^ 
v = N mod 2. When lattice effects 2 ^ are taken into 
account, inter- valley exchange weakly favors the layer- 
antifcrromagnet (LAF)£ state. (That observation does 
not definitively make the case for this particular state, 
since inter-valley exchange could be less important than 
correlation effects.) For a LAF state, cth = even though 
time-reversal symmetry is broken by opposite spin polar- 
izations on the top and bottom layers. In a continuum 
model that neglects lattice effects and has SU(4) sym- 
metry, the three states are all members of the continuous 
family of degenerate broken symmetry states that is sig- 
naled by divergence of x| z . 

In addition to their difference in edge state proper- 
ties, the QAH, QSH and LAF states also exhibit different 
responses 2 ^ to Zeeman fields which couple to spin and can 
be realized in the absence of orbital coupling by applying 
magnetic fields parallel to graphene layers. In the QAH 
case, the ground state is unchanged but the quasiparti- 
clc gap is reduced, vanishing when the Zeeman-coupling 
strength is equal to the ground state gap via a mecha- 
nism reminiscent of the Clogston limit in superconduc- 
tors. The QSH and LAF states respond to Zeeman fields 
by establishing noncollinear spin states within each valley 
and evolving toward an unusual kind of exciton conden- 
sate in the strong Zeeman-coupling limit. The quasipar- 
ticle gap of QSH and LAF states is however independent 
of Zeeman-coupling strength drawing a sharp distinction 
with the QAH case. It is therefore possible to use Zee- 
man responses and edge state signatures to identify the 
character of the bilayer ground state experimentally^ 
On the other hand, the three states respond very simi- 
larly to an electric field between the layers, which can 
induce first order transitions at which the total layer 
polarization jumpsJ^^ This state favored by an inter- 



layer electric field is a quantum valley Hall state with 
zero Hall conductivity but, in the absence of scattering 
which mixes valleys counterpropagating valley-resolved 
edge states. We note that a small electric field between 
the layers can also possibly stabilize a state in which one 
flavor is polarized in a sense opposite to the other three 
and charge, valley, and spin Hall conductivities are all 
nonzero £ 

One consequence of the above mentioned spontaneous 
quantum Hall effects^ in each valley is that the gapped 
states have a very simple adiabatic evolution with mag- 
netic field B in which the gap remains open, but the car- 
rier density at which the gap appears varies with field. 
The total filling factor inside the gap at finite B is v 
for a state with a Hall conductance equal to ve 2 /h at 
zero field. The topological character— of the gapped 
states in BLG is induced^ by weak bare electron-electron 
interactions and is enhanced by the nearly flat bands 
with chiral pseudospin textures, and is distinct from the 
Z2 topological insulators^ 2 - - — that require strong intrin- 
sic spin-orbital coupling. Although interaction-induced 
topological states apparently do not occur in single-layer 
graphene, as proposed^ on the basis of extended Hub- 
bard models with at least finite interaction strength, it 
appears quite possible that they do occur in bilayers via 
a weak interaction instability^ 

There are three previous article a 8 ' 17 ' 22 which discuss 
similar RG calculations. Ref. [l?] reports equivalent RG 
flow equations for the SU (4) invariant model in terms 
of interaction parameters that are related to ours by 
9o.-z = Ts ± Td and 2g± = Fx ■ Several references have 
started their RG flows from effective interactions which 
include symmetry-allowed corrections due to lattice ef- 
fects. Ref. [22] attempts to account for the long-range 
of the Coulomb interaction by explicitly accounting for 
screening and then replacing the Coulomb amplitude by 
an average value, and estimates the ground state by ap- 
plying mean-field theory to interactions from PRG flows 
that have been truncated at an intermediate interaction 
strength. All calculations recognize the competition be- 
tween nematic and gapped states, but small differences 
in approximations or in estimates of bare coupling pa- 
rameters have lead to different final conclusions. We 
note that the recent elaboratio n] 20 - 34 of the calculations 
which motivated nematic state proposal o 17 ' 22 recognize 
that gapped states^ are also a possibility. 

Although these calculations, like ours, shed light on 
the physics at play in determining the ground state, 
it is difficult to make definitive conclusions based on 
theory alone and we must ultimately appeal to exper- 
iment. Here the situation is also confusing. Some ex- 
perimental studies 1 ^ - — of high-quality suspended bilay- 
ers appear to reveal gaps whose size, carrier-density- 
dependence, and temperature-dependence, is roughly 
consistent with mean-field-theory predictions . 14 ' 36 These 
observations are consistent with the conclusions of this 
paper. On the other hand this behavior is not observed 
in all samples, even in samples which appear to be sim- 
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ilar to the gapped samples. One experimental study- 1 - 1 - 
focuses on the magnetic-field B dependence of the gap 
at filling factor v = 4 in a sample which does not show a 
gap at B = 0. The anomalous weak B persistence of the 
v = 4 gap in this experiment appears to be consistent 
with earlier studies by Yacoby and collaborators.^ 1 ^ Al- 
though interpreted as evidence for the nematic study, the 
persistence of the v — 4 gap is, it appears to us, equally 
consistent with a field-driven crossover or transition to a 
gapped oh = 4e 2 //i state. Although this gap is likely not 
the ground state at B = it is favored at filling factor 
v = 4, because its gap is then pinned to the Fermi level. 
In this interpretation the increased conductance observed 
at weak fields in Ref. [n] would be associated with a first 
order transitio n 14 ! 35 between ch = 4e 2 //i and cth = 
states which induces a network of current-supporting 5 ^ 
domain walls.— 

Similar spontaneously broken symmetry physics^ is 
also likely in thicker high quality graphene films with chi- 
ral (ABC) stacking order J±££ In an A^-layer chiral stackj 4 - 
the low-energy sublattice sites are localized in the outer- 
most layers, A\ and B^. Hopping occurs between these 
two sites via an TV-step virtual process which leads to p N 
energy dispersions and an Nn Berry phase. For larger N, 
the low-energy bands are increasingly flat and the pseu- 



dospin chirality larger, at least when weak remote hop- 
ping processes are neglected, leading to larger opening 
for many-body interaction effects. In the simplified chiral 
model, the density of states D(E) ~ ij( 2 ~ N )/ N diverges 
as E approaches zero for N > 2 whereas it remains finite 
for N = 2. In the PRG language, this difference corre- 
sponds to the fact that the short-range interactions at 
tree level are marginal for N — 2 but relevant for TV > 2. 
Since Vx — at tree level, as discussed in Section III-C, 
it seems that gapped spontaneous quantum Hall states^ 
which break layer-inversion symmetry and produce large 
momentum space Berry curvatures can occur in multi- 
layers as well. Future experimental and theoretical work 
will be necessary to sort out the competition between in- 
teractions and remote hopping terms in the Hamiltonian. 
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